//visRayCasting.cpp

#include "visRayCasting.h"
namespace visualization
{
void visRayCasting::setImageSize(int size)
{
	if(size < 20)
	{
		cout << "AT LINE " << __LINE__ << " OF FILE " << __FILE__ << endl;
		cout << "Seems that the parameter for the size of the image is too small(" << size << ")" << endl;
		exit(-1);
	}
	imgSize = size;
	
	rayCastingImage = new image*[imgSize];
	if(!rayCastingImage)
	{
		cout << "AT LINE " << __LINE__ << " OF FILE " << __FILE__ << endl;
		cout << "No memory for it!" << endl;
		exit(-1);
	}
	for(int i = 0; i < imgSize; i++)
	{
		rayCastingImage[i] = new image[imgSize];
		if(!rayCastingImage[i])
		{
			cout << "AT LINE " << __LINE__ << " OF FILE " << __FILE__ << endl;
			cout << "No memory for it!" << endl;
			exit(-1);
		}
	}
}

void visRayCasting::computeBoundingCube()
{
	double xmin, xmax, ymin, ymax, zmin, zmax;
	dataField -> getBound(xmin, xmax, ymin, ymax, zmin, zmax);

	cout << "xmax: " << xmax << endl;
	cout << "xmin: " << xmin << endl;

	cout << "ymax: " << ymax << endl;
	cout << "ymin: " << ymin << endl;
	
	cout << "zmax: " << zmax << endl;
	cout << "zmin: " << zmin << endl;

	// The center of the cube.
	center[0] = (xmax + xmin)/2;
	center[1] = (ymax + ymin)/2;
	center[2] = (zmax + zmin)/2;

	double xR = xmax - xmin;
	double yR = ymax - ymin;
	double zR = zmax - zmin;

	// Get the edge length of the cube. Note that the edge_length is the diameter of the 
	// circumcircle of the cube surrounding the data field.
	edge_length = sqrt(xR * xR + yR * yR + zR * zR);

	cout << "center: " << center[0] << "	" << center[1] << "	" << center[2] << endl;
	cout << "edge_len: " << edge_length << endl;

}

void visRayCasting::setViewingPosition(double position[], int dim)
{
	// Compute the center and edge_length here.
	computeBoundingCube();

	if(dim != 3)
	{
		cout << "The parameter for the light direction should have three components! " << endl;
		exit(-1);
	}
	// In the current design,
	// If it coincides with the center of the data field, then it will be impossible to
	if( abs(center[0]-position[0])<1e-6 && abs(center[1]-position[1])<1e-6 && abs(center[1]-position[2])<1e-6)
	{
		cout << "The viewing position coincides with the center of the center of the data field!"
			 << "Change another direction parameter please!" << endl;
		exit(-1);
	}
	for(int i = 0; i < dim; i++)
	{
		direction_to_origin[i] = position[i];
	}
}

void visRayCasting::computeDirectionToCenter()
{
	direction_to_center[0] = center[0] - direction_to_origin[0];
	direction_to_center[1] = center[1] - direction_to_origin[1];
	direction_to_center[2] = center[2] - direction_to_origin[2];
}

void visRayCasting::resample()
{
	// 1. Compute the plane(in fact a square) that is orthogonal to the light direction,
	//    which points to the center of the cube. The result is four points.
	// 2. According to the image size, it can be determined where the light enters the
	//    data field.
	// 3. According to the "samplingRate", interpolate to get the value at the correspo-
	//    nding resample points.
	
	// 1
	// The direciton "direction_to_center" will be needed to calculate the four points.
	computeDirectionToCenter();
	
	// Normalize it !
	normalize(direction_to_center, direction_to_center);

	// calculate an arbitrary direction that is orthogonal to "direction_to_center".
	double tmpDir1[3], tmpDir2[3];
	getOrthogonalDirection(direction_to_center, tmpDir1);
	normalize(tmpDir1, tmpDir1);

	crossMultiply(direction_to_center, tmpDir1, tmpDir2);
	normalize(tmpDir2, tmpDir2);

	/*  1         2
	 *  *---------*
	 *  |         |
	 *  |    *c   |
	 *  |         |
	 *  *---------*
	 *  3         4
	*/

	// Coordinates of the four vertexes.
	double v1[3], v2[3], v3[3], v4[3];
	double diag = edge_length/sqrt(2.0);

	advanceAPoint(center, tmpDir1,  1, diag, v1);
	advanceAPoint(center, tmpDir1, -1, diag, v4);
	advanceAPoint(center, tmpDir2,  1, diag, v2);
	advanceAPoint(center, tmpDir2, -1, diag, v3);

	// Now the coordinates the vertexes are got.
	advanceAPoint(v1, direction_to_center, -1, edge_length/2, v1);
	advanceAPoint(v2, direction_to_center, -1, edge_length/2, v2);
	advanceAPoint(v3, direction_to_center, -1, edge_length/2, v3);
	advanceAPoint(v4, direction_to_center, -1, edge_length/2, v4);

	cout << "direction_to_center: " << direction_to_center[0] << " " << direction_to_center[1] << " " << direction_to_center[2] << endl;

	// 2
	// The entrances of the light can be determined now by the image size "imageSize".
	// Then for each of them, the resampling can be carried on.
	double dirV1ToV2[3] = {
							 tmpDir2[0] - tmpDir1[0],
							 tmpDir2[1] - tmpDir1[1],
							 tmpDir2[2] - tmpDir1[2]
						  };
	normalize(dirV1ToV2, dirV1ToV2);
	double dirV1ToV3[3] = {
							 -tmpDir2[0] - tmpDir1[0],
							 -tmpDir2[1] - tmpDir1[1],
							 -tmpDir2[2] - tmpDir1[2]
						  };
	normalize(dirV1ToV3, dirV1ToV3);
/*		for(each row)
		for(each col)
		{
			1. Resample according to samplingRate to get sample point with location.
			2. Invoke resample and combine function to 
			   Advance the point along the eye light to go through the data field. 
			   1. For each sample point:
				  1. interpolate to get value of the sample point.
				  2. according to the value, get corresponding color and opacity.
			   2. Combine the results to make a pixel.
		}
*/
	// This is the step length, that is, the density of the image.		
	double density = edge_length / (imgSize - 1);

	// The interval between two neighbouring resampling points.
	double samplingInterval = edge_length / (smpRate - 1);

	// The point at the beginning of each row.
	double tmpP[3];

	// If later the sampling rate is gona to change.
	// The sample points are saved in pntsOnLight one ray after another. The parameter 3 is the dimension 
	// of a point.
	const int dimension = 3;
	const int pointN = smpRate*imgSize*imgSize;
	pntsOnLight = new double*[pointN];
	if(!pntsOnLight)
	{
		cout << "AT LINE " << __LINE__ << " OF FILE " << __FILE__ << endl;
		cout << "SORRY,NO MEMORYR! " << endl;
		exit(-1);
	}
	for(int i = 0; i < pointN; i++)
	{
		pntsOnLight[i] = new double[dimension];
		if(!pntsOnLight[i])
		{
			cout << "AT LINE " << __LINE__ << " OF FILE " << __FILE__ << endl;
			cout << "SORRY, NO MEMORY!" << endl;
			exit(-1);
		}
	}
	sampledPoints = new double[pointN];
	if(!sampledPoints)
	{
		cout << "AT LINE " <<  __LINE__ << " OF FILE " << __FILE__ << endl;
		cout << "SORRY, NO MOMERY!" << endl;
		exit(-1);
	}

	int pN = 0;
	// Compute the positions of sampling points.
	for(int v1ToV3 = 0; v1ToV3 < imgSize; v1ToV3++)
	{
		advanceAPoint(v1, dirV1ToV3, 1, v1ToV3*density, tmpP);

		for(int v1ToV2 = 0; v1ToV2 < imgSize; v1ToV2++)
		{
			// The wanted initial location.
			advanceAPoint(tmpP, dirV1ToV2, 1, v1ToV2*density, pntsOnLight[pN]);

			// Advance the point along the eye light to get resamplingRate resampling points.
			for(int i = 1; i < smpRate; i++)
			{
				advanceAPoint(pntsOnLight[pN], direction_to_center, 1, samplingInterval, pntsOnLight[pN+1]);
				pN++;
			}
			// This is important. This line of code has cost almost one day.
			pN++;
		}
	}
	// Interpolate to get the value at the position of each sampling point.
	std::vector<double> value;
	for(int i = 0; i < pointN; i++)
	{
		value.clear();
		if( !(dataField -> interpolatedValue(pntsOnLight[i], 3, &value, 5, 5)) )
		{
			// If it can not be interpolated, the sampledPoints would be the maximum of the double type.
			sampledPoints[i] = DBL_MAX;
			continue;
		}
		sampledPoints[i] = value[0];
	}
	cout << "Interpolation ok!" << endl;
}

void visRayCasting::getOrthogonalDirection(double direction[], double tmpDir[])
{
	char x = abs(direction[0]) < 1e-10 ? 0 : 1;
	char y = abs(direction[1]) < 1e-10 ? 0 : 1;
	char z = abs(direction[2]) < 1e-10 ? 0 : 1;
	int flag = (x << 2) + (y << 1) + z;

	switch(flag)
	{
		case 0:
		{
			// (x, y, z) = (0, 0, 0)
			cout << "AT LINE " << __LINE__ << " OF FILE " << __FILE__ << endl;
			cout << "All components are zeros, so the orthogonal direction can't be calculated."
				 << endl;
			exit(-1);
		}
		case 7:
		{
			// (x, y, z) = (1, 1, 1)
			tmpDir[0] = tmpDir[1] = 1;
			tmpDir[2] = -(direction[0] + direction[1]) / direction[2];
			break;
		}

		case 1:
		{
			// (x, y, z) = (0, 0, 1)
			tmpDir[0] = tmpDir[1] = 1; 
			tmpDir[2] = 0;
			break;
		}
		case 2:
		{
			// (x, y, z) = (0, 1, 0)
			tmpDir[0] = tmpDir[2] = 1;
			tmpDir[1] = 0;
			break;
		}
		case 4:
		{
			// (x, y, z) = (1, 0, 0)
			tmpDir[1] = tmpDir[2] = 1;
			tmpDir[0] = 0;
			break;
		}

		case 3:
		{
			// (x, y, z) = (0, 1, 1)
			tmpDir[0] = 0;
			tmpDir[1] = 1;
			tmpDir[2] = -direction[1]/direction[2];
			break;
		}
		case 5:
		{
			// (x, y, z) = (1, 0, 1)
			tmpDir[0] = 1;
			tmpDir[1] = 0;
			tmpDir[2] = -direction[0]/direction[2];
			break;
		}
		case 6:
		{
			// (x, y, z) = (1, 1, 0)
			tmpDir[0] = 1;
			tmpDir[1] = -direction[0]/direction[1];
			tmpDir[2] = 0;
			break;
		}
	}// switch
	normalize(tmpDir, tmpDir);
}// getOrthogonalDirection

void visRayCasting::crossMultiply(double direction1[], double direction2[], double result[])
{
	// x1: direction1[0]
	// y1: direction1[1]
	// z1: direction1[2]

	// x2: direction2[0]
	// y2: direction2[1]
	// z2: direction2[2]

	// y1 * z2 - y2 * z1
	result[0] = direction1[1] * direction2[2] - direction2[1] * direction1[2];
	// x2 * z1 - x1 * z2
	result[1] = direction2[0] * direction1[2] - direction1[0] * direction2[2];
	// x1 * y2 - x2 * y1
	result[2] = direction1[0] * direction2[1] - direction2[0] * direction1[1];

	normalize(result, result);
}

void visRayCasting::normalize(double vec[3], double result[3])
{
	double sum = pow(vec[0], 2) + pow(vec[1], 2) + pow(vec[2], 2);
	sum = sqrt(sum);
	if(abs(sum) < 1e-10)
	{
		cout << "AT LINE " << __LINE__ << " OF FILE " << __FILE__ << endl;
		cout << "Program stops because some number is going to be devided by 0. " << endl;
		exit(-1);
	}
	for(int i = 0; i < 3; i++)
	{
		result[i] = vec[i] / sum;
	}
}

inline void visRayCasting::advanceAPoint(const double point[], const double direction[], const double flag, const double length, double result[3])
{
	for(int i = 0; i < 3; i++)
	{
		result[i] = point[i] + direction[i] * length * flag;
	}
}
bool visRayCasting::computeGradient(double &component, const double step, const int flag, const double centralPoint[3])
{
	double forwardP[3], backwardP[3];
	for(int i = 0; i < 3; i++)
	{
		forwardP[i] = backwardP[i] = centralPoint[i];
	}
	std::vector<double> fvalue, bvalue;
	int index;

	switch(flag)
	{
		// X
		case 0:
		{
			index = 0;
			break;
		}
		// Y
		case 1:
		{
			index = 1;
			break;
		}
		// Z
		case 2:
		{
			index = 2;
			break;
		}
	}
	forwardP[index] = centralPoint[index] + step;
	backwardP[index] = centralPoint[index] - step;

	if( !dataField->interpolatedValue(forwardP, 3, &fvalue, 10, 10) 
		|| !dataField->interpolatedValue(backwardP, 3, &bvalue, 10, 10) )
	{
		return 0;
	}

	component = (fvalue[0] - bvalue[0])/(2*step);

	return 1;
}

void visRayCasting::combineToFormImage(const int direction)
{
	cout << "Begin to combine image.." << endl;
	//initialize the image.
	for(int row = 0; row < imgSize; row++)
	{
		for(int col = 0; col < imgSize; col++)
		{
			for(int i = 0; i < 3; i++)
			{
				rayCastingImage[row][col][i] = i * 60 + 100;
			}
		}
	}

	// Combine in a way from back to front.
	if(1 == direction)
	{
		int j = 0;
		for(int row = 0; row < imgSize; row++)
		{
			for(int col = 0; col < imgSize; col++)
			{
				for(int i = j + smpRate - 1; i >= j; i--)
				{
					if(sampledPoints[i] < DBL_MAX/10.)
					{
						double *colorAndOpacity = transFunc.getColorAndOpacity(sampledPoints[i]);
						// For each color component(RGB).
						for(int colorCom = 0; colorCom < 3; colorCom++)
						{
							rayCastingImage[row][col][colorCom] = (unsigned char)( rayCastingImage[row][col][colorCom] * ( 1 - colorAndOpacity[3] ) +
															   	  colorAndOpacity[colorCom] * colorAndOpacity[3] );
						}
					}
				}
				j += smpRate;
			}
		}
	}
	// Combine in a way from front to back.
	else if(0 == direction)
	{
		for(int row = 0; row < imgSize; row++)
		{
			for(int col = 0; col < imgSize; col++)
			{
				for(int i = 0; i < 3; i++)
				{
					rayCastingImage[row][col][i] = 0;
				}
			}
		}

		double t = 1;
		int j = 0;
		double transparacyThreshold = 0.01;

int at = 150;

		for(int row = 0; row < imgSize; row++)
		{
			for(int col = 0; col < imgSize; col++)
			{
				int peel = 1;
				for(int i = j; (i <= j + smpRate - 1) /*&& (t > transparacyThreshold)*/; i++)
				{
					if(t < transparacyThreshold)
					{
						if(at > peel)
						{
							t = 1;
							peel++;
						}
						else
						{
							break;
						}
					}
					if(sampledPoints[i] < DBL_MAX/10.)
					{
						double *colorAndOpacity = transFunc.getColorAndOpacity(sampledPoints[i]);
						// For each color component(RGB).
if(at == peel)
{						
						for(int colorCom = 0; colorCom < 3; colorCom++)
						{
							rayCastingImage[row][col][colorCom] += (unsigned char)(t * colorAndOpacity[colorCom]);
						}
}
						t *= (1 - colorAndOpacity[3]);
					}
				} // for
				j += smpRate;
			}
		}

		for(int row = 0; row < imgSize; row++)
		{
			for(int col = 0; col < imgSize; col++)
			{
				for(int i = 0; i < 3; i++)
				{
					rayCastingImage[row][col][i] += (unsigned char)(t * (i * 60));
				}
			}
		}
	}
	// This function explores the first derivative property...
/*	else if(-1 == dir)
	{
		// 1. For each point, compute the first derivative, or gradient.
		// 2. According to the magnitude of the gradient, assign color and opacity for 
		//    each point.
		for(int i = 0; i < 3; i++)			
		{
			generatedImg[row][col][i] = i * 50 + 60;
		}
		
		// Three components of the gradient.
		double gx, gy, gz;
		double magnitude;
		double step = 0.5;

*/
/*		double forward[3], backward[3];

		std::vector<double> valueF, valueB;
*/
/*		for(int i = 0; i < smpRate; i++)
		{
*/

			/*if( !dataField -> interpolatedValue(pntOnLight[i], 3, &value, 10, 10)) )
			{
				continue;
			}
			*/


/*			// For each point, get gx, gy and gz, and then assign color and opacity to the point.
			if( !computeGradient(gx, step, 0, pntOnLight[i])
			   || !computeGradient(gy, step, 1, pntOnLight[i])
			   || !computeGradient(gz, step, 2, pntOnLight[i])
			  )
			{
				continue;
			}
			magnitude = sqrt( pow(gx, 2) + pow(gy, 2) + pow(gz, 2) );

			double *colorAndOpacity = transFunc.getColorAndOpacity(magnitude);

			// For each color component(RGB).
			for(int colorCom = 0; colorCom < 3; colorCom++)
			{
				generatedImg[row][col][colorCom] = generatedImg[row][col][colorCom] * ( 1 - colorAndOpacity[3] ) +
												   colorAndOpacity[colorCom] * colorAndOpacity[3];
			}
		}
 		
	}
*/	else
	{
		cout << "AT LINE " << __LINE__ << " OF FILE " << __FILE__ << endl;
		cout << "Correct combination way should be selected! " << endl;
		exit(-1);
	}
}

}

